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ABSTRACT 

This paper extends in two directions the results of prior work on generalized linear covariance analysis of 
both batch least-squares and sequential estimators. The first is an improved treatment of process noise in the 
batch, or epoch state, estimator with an epoch time that may be later than some or all of the measurements 
in the batch. The second is to account for process noise in specifying the gains in the epoch state estimator. 
We establish the conditions under which the latter estimator is equivalent to the Kalman filter. 

INTRODUCTION 

Reference 1 gave an overview of a comprehensive approach to 1 i near covariance analysis, closely follow- 
ing the approach of Refs. 2 and 3, and including references to earlier studies of this subject. The models 
presented include both batch least-squares methods, which originated with Gauss, 4 and sequential estima- 
tors, which were pioneered by Kalman. 5 Both batch least-squares estimation and Kalman filtering methods 
produce the covariance matrix encapsulating both the expected errors of the estimated parameters and the 
correlations among them. The number of parameters describing a dynamical system can be quite large, 
and it is often impractical to solve for all of them, especially in real-time applications. “Consider covari- 
ance analysis,” the subject of Refs. 1-3, treats estimators that solve for only a subset of the parameters, 
but considers the effect on the covariance of the parameters that are not solved for. Our approach differs 
from the Schmidt Kalman filter, 6 which takes the consider parameters into account in computing the up- 
date gains, but updates neither their values nor their covariances. Our consider parameters, in contrast, arc 
completely ignored by the estimation algorithm. Their effect on the accuracy of the estimates obtained is 
considered in the covariance analysis, though, accounting for their name. This treatment is very similar to 
that described in Bierman 7 and in Tapley, et al. 8 The method described here follows Refs. 1-3 and differs 
from many other treatments in explicitly partitioning the errors into subspaces containing only the influence 
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of the measurement noise, process noise, and a priori solve-for and consider covariances. McReynolds 9 
has also implemented such a partitioning in both batch and sequential filter/smoother error budget analysis. 
Explicitly exhibiting the a priori component of the covariance is useful for observability analysis because a 
priori information about states that arc observable is “forgotten” as more data arc processed. 10 

The treatment of the batch estimation in Refs. 1-3 assumed that the estimator did not account for process 
noise in computing the update gains, as is customary in batch least-squares estimation. The estimators 
considered in this paper, which we call epoch-state estimators to distinguish them from the usual batch 
estimators, take process noise into consideration. We show that the usual batch estimator is recovered when 
process noise is ignored by the estimator. Process noise may or may not be included in the estimation 
process, but it is always accounted for in the linear covariance analysis, as in Refs. 1-3, because it can 
contribute to estimation errors whether the estimator knows about it or not. 

We define the epoch time to be the time at which the a priori information is specified and at which the state 
update is performed. The analysis of the batch estimator in Refs. 2 and 3 implicitly assumed that the epoch 
time is prior to all of the measurements, which is the usual case, and which the term a priori implies. In this 
case the a priori errors at the update time arc uncorrelated with the process noise during the measurement 
span. It is sometimes convenient to allow the epoch time to be later than some or all of the measurements, 
though, and Ref. 1 considered this case. It continued to assume zero correlations between the a priori errors 
and the process noise, however. It may be more reasonable to assume that the a priori errors include the 
process noise accumulated between the beginning of the observation span and the epoch time, if the latter is 
later than some or all of the measurements. This modifies the manner in which process noise appeal's in the 
covariance analysis of the epoch-state estimator, and is the principal subject of this paper. 

To make this paper self-contained and to establish notation, the first several sections of this paper contain 
a review of the material in Ref. 1. We have, however, attempted to eliminate unnecessary duplication with 
the earlier paper. We first consider the linear state dynamics models, and the partitioning of state errors into 
those arising from a priori errors, from measurement errors, and from process noise. The dynamics does not 
distinguish between solve-for and consider parameters, so this distinction is not made until we turn to the 
processing of measurement information. The measurement update equations involve the projection matrix 
onto the solve-for subspace, because the update does not affect the consider parameters. We then compute 
the contributions from a priori errors, measurement errors, and process noise to the covariance matrix of 
the full state. The measurement errors are uncorrelated with both the a priori errors and the process noise, 
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so the contribution of the measurement errors to the covariance can be defined unambiguously. However, 
we now account for correlations between the a priori errors and the process noise, and we include these 
correlations in the process noise paid of the covariance. 

We next discuss the selection of the update gains to minimize the covariance of the estimation errors at 
the epoch time, taking into account a model of the process noise in the gain specification. We show that 
our results recover the usual batch least-squares filter when the process noise is zero. If process noise is 
not neglected, the epoch-state estimators fall into three classes: 1) when no measurements follow the epoch 
time, the estimator can be operated sequentially, but the covariance analysis of Refs. 1-3 must be modified; 
2) when no measurements precede the epoch time, the covariance analysis of Refs. 1-3 is valid, but the 
estimator cannot be operated sequentially; 3) when some measurements precede and some follow the epoch 
time, sequential operation of the estimator is precluded, and the covariance analysis of Refs. 1-3 must be 
modified. 

The sequential formulation of this epoch-state estimator has been applied to an analysis of satellite close 
approaches. 11 For this application, two filters are run concurrently, one constrained by the condition that the 
distance of closest approach is greater than some effective hard body radius, and the other constrained so 
that the closest approach distance is less than the hard body radius. The constraint is implemented using the 
technique of Zanetti, et al. 12 An epoch state filter with epoch time in the future is required because the time 
of closest approach, where the constrained update is enforced, is in the future of all the measurements. 

We close the paper by showing that the sequential epoch-state estimator is equivalent to the Kalman filter. 5 
and by presenting a summary and our conclusions. 

MODELING APPROACH 
General Dynamics Model 

The state vector x is an //-dimensional vector of parameters that completely characterizes a system, in- 
cluding both time-dependent and constant parameters. By assumption, the state vector evolves according 
to 

^x(t) = f(x(t),t) + w(t) (1) 

where the dynamic noise w(t) is a Gaussian white noise process with mean and covariance given by 

E[ w (f)]=O n and E [w(i)w(r) T ] = Q(t)5(t - r), (2) 
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with E [...] denoting the expectation operator and 0 n denoting an n-dimensional vector of zeros. In the 
covariance equation, Q is the n x n dynamic noise spectral density matrix and 5{t — t) denotes the Dirac 
delta, or unit impulse, function. 


The true value of the state vector is never exactly known, but can only be estimated. The state estimate 


vector, x, evolves according to 


The state error vector, given by 


-±{t) = f(x(t),t). 


e(t) = x(t) -x(t), 


is assumed to always remain small. Then first-order linear error analysis techniques, as commonly found in 


an extended Kalman filter, give 


—e(t) = A(t)e(t) + w(i), 


where 


A(t) = 


df (x(t),t) 
dx(t) 


Formal integration of equation (5) gives 


e(t) = t')e(t') + w d (t, t') 


where e(f') the eiTor vector at time t', the state transition matrix, <b(i, t'), is the solution of 


= AmM 


with the initial condition 


t 1 ) = /„ = the n x n identity matrix, 


and the random excitation vector, w d(t, f), is given by the formal integral 


w d(t,t')= [ $(f,r)w(r)dr. 
Jt' 


We will never need to evaluate this mathematically ill-defined integral, but we use it to formally compute 
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covariances, which are the well-defined quantities of interest. 

The state vector error at time t can be written as the sum of three contributions: the error e a (t), due to an 
a priori error at an epoch time f*. the error e v (t) due to noise in the measurements used for the estimation, 
and the error e w (t) due to dynamic noise. 

e(f) = e a (t) + e v (t) + e w (t). (11) 

Equation (7) can be used to propagate the components of the state error vector from the epoch time f* to 
some other time t: 


e a (t) = <f*(f,t*)e a (f*), 

(12a) 

e v (t) = $(t,t*)e v (t*), 

(12b) 

e w (t) = 4>(t, t*)e w (t*) + w d (t, U). 

(12c) 


General Estimation Model 

It is usually not necessary to produce estimates of all of the state parameters, but only of a subset of 
n s solve-for parameters. The remaining n c = n — n s parameters arc called consider parameters because 
they arc not estimated but contain uncertainties that must be considered in the error analysis. The solve-for 
parameter vector s (t) and the consider parameter vector c (t) arc given by 

s (t) = S(t)x.(t) and c (t) = C(t)x(f) (13) 


where the n s x n matrix S (t) and the n c x n matrix C it) arc such that the matrix 


M = 


S 


C 


(14) 


is non-singular*. The inverse of M is partitioned into an n x n s matrix S and an n x n c matrix C: 


M~ 1 


S, C . 


( 15 ) 


*Here and subsequently we suppress the time arguments where there is no ambiguity so as to simplify the notation. 
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The properties of the matrix inverse then lead immediately to the identities 


SS = I ns , CC = I Uc , SC = 0 naXnc , CS = O ncXns , (16) 

and 

SS + CC = I n . (17) 

In the usual case that the elements of the solve-for and consider vectors are merely selected and possibly 
permuted components of the state vector, the matrix M is an orthogonal permutation matrix. In this case, 
and in any case for which M is orthogonal, the matrices S and C arc just the transposes of S and C, 
respectively, which makes inversion of M unnecessary and simplifies many of the following equations. 

It follows from equations (13) and (17) that 

x(t) = S(t)s(t) + C(t)c(t). (18) 

We also have an estimated value s(t) of the solve-for vector and an assumed value c(t) of consider vector. 
Both are propagated by the dynamics, but only the solve-for vector is updated by the estimator. Relations 
similar to equation (13) give the estimated solve-for vector and the assumed consider vector the in terms of 
the estimate of the full state x(t). Thus, the error vector can also be expressed in terms of its solve-for and 
consider components as 


e s (t) = s (t) — s(t) = S(t)e(t ) (19a) 

e c (t) = c (t) — c(t) = C(t)e(t) (19b) 

and 

e(f) = S(t)e s (t) + C(t)e c (t). (20) 

The subscripts s and c arc independent of the subscripts a, v, and w in the previous section. Both the solve- 
for errors and the consider errors have contributions arising from the a priori estimate, the measurement 
noise, and the process noise. We will avoid double subscripts, however. 

An estimator produces its state estimates based on an a priori estimate at an epoch time t* and infor- 
mation obtained from measurements obtained at discrete times t t . Let y, be an mi-dimensional vector of 
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measurement values obtained at times t t . Measurements are related to the state vector by the following 
measurement model: 

Yi = hj(xj,fj) + Vj, (21) 

where v* is a Gaussian white noise process with mean and covariance given by 

E [vj] = 0 , E [vjvf ] = Ri, and E [v, t vj] = 0, for i / j. (22) 

The functions h ; are assumed to be known functions of imprecisely known arguments. Therefore, it is 
possible to compute predicted measurements from 


y* = h *(x ? : ,U) (23) 

where x“ is the state estimate before incorporating the measurement information. Conversely, x- h denotes 
the state estimate at time t{ including the measurement information.^ A Taylor series expansion gives the 
difference between the observed and computed measurements, or “o-minus-c,” to first order in e, as 

r- = y i- y i = Hie~ + Vjj (24) 

where the measurement sensitivity matrix Hi is given by 

dhi(xi,ti) 

Ilj = — 

9x, : 

This linearization is familial - from the extended Kalman filter. The o-minus-c before the update, r“, is 
commonly referred to as the “innovation,” and the o-minus-c after the update as the “residual,” though this 
terminology is by no means universal. 


(25) 


^We use this notation in the sequel to denote a priori and a posteriori values of various quantities. 
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Epoch State Estimator 


An epoch-state estimator produces an estimate by performing a single update at the epoch time f* based 
on a single batch measurement y that combines the N measurements y, made at various times. Thus, 


yi 


yi 


r i 


A 


A 


y n 

, y = 

y n 

li 

>> 

i 

II 



(26) 


The update is given by: 


N 

S+ = S" + K r = s~ + ^2 


(27) 


i — 1 

where s“ is the a priori estimate of s(t*) and K is a gain matrix consisting of a “row” of gain matrices for 
each measurement, /ij. Substituting equations (24) and (7), with t = t, and t! = f*, into equation (26) gives 


r = + u d + v 


(28) 


where 





1 

* 

■+o 

-+o 

1 


Ufil 


H 1 vr d (ti, f*) 


Vl 

H = 

Hn,* 


77/v<b(t/v, t*) 

, u d = 

UdN 


H N w d (t N ,t *) 

, v = 

VAT 


(29) 


The state estimates must be propagated from the epoch time t* to the measurement times t, in order to 
compute the measurement sensitivity matrices H, by equation (25), and these arc propagated back to the 
epoch time to obtain the matrices H,„. 


Since the estimation process does not alter the consider parameters, the estimation error after an update 
is related to its value before the update by 

N 

= e“ + 5* (e+ - e“ ) = e“ - 5* ^ K t (H^e~ + v* + u di ) , (30) 

i = 1 


where S< = S(t *). Now equation (11) is used to partition the state error into a priori, measurement, and 



process noise components: 


e 


+ 

a* 


e 


+ 

V* 


e 


+ 

w* 


N 


I n -S*J2 K iH i: 


i= 1 


N 

-S'* ^ K,y t 
2—1 
N 

A ^ " KiUdi- 
i = 1 


(31a) 

(31b) 

(31c) 


COVARIANCE ANALYSIS 

The function of an estimator is to produce an estimate s+ given measurements y,. Error analysis, however, 
does not actually compute an estimate, but determines how good an estimate would be produced in a given 
situation. The covariance matrix of the estimation error, defined, assuming E [e] = 0, by 


P(t*) = E [e(t*)e(f*) T ] , (32) 

provides a statistical measure of the quality of an estimate at the epoch time t* of a given scenario^. The a 
priori covariance is simply 

P~= E[e-(e-) T ] (33) 

The post-update ( a posteriori) covariance is obtained by substituting equation (31) into 

P+ = E [e+(e+) T ] = E [(e+ + e+ + e+J(e+ + e+ + e+J T ] (34) 

This computation is complicated by correlations between e+, and e^*. To see how these arise, consider a 
time to (not the epoch time) prior to the epoch time t* and to all of the measurement times t t . Equation (7) 
relates the a priori error at the epoch time f* to the error at to by 


e* = $(t*,f 0 )e 0 + w d (U,t 0 ). (35) 

Reference 1 assumed that c) is not coiTelated with e^, but it is more natural to assume that eg is not 

correlated with with the result that the second term on the right side of equation (35) will give rise 

Tt should be kept in mind that if systematic errors cause E [e] ^ 0. then E(t*) will describe the mean square enor matrix, 
rather than the covariance. 
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to correlations between e“ and e+* unless the epoch time is earlier than all the measurement tiniest This 
is where the “arrow of time” enters the picture. The a priori error c“ depends on the process noise in the 
past of t#, but not on process noise in the future. We will see that the correlations introduced only involve 
integrals of the process noise power spectral density over the overlaps of the interval between i* and to and 
the intervals between t* and ti, so the time to will not appear in our final equations. 

The covariance matrix can be written as the sum of a priori, measurement noise, and dynamic noise 
contributions. We include the correlation between e+ and e+* in the process noise part, giving 


P t = P a* + P v* + P w * > where 

(36a) 

P a* = E [e+ (e+ ) T ] , 

(36b) 

P+ = E [e+ (e+ ) T ] , 

(36c) 

p tt* = E [e+*(e+J T ] + E [e+,(e+ ) T ] + E [e+ (e+*) T ] . 

(36d) 


The a priori error is uncorrelated with the measurement errors because the a priori estimate is assumed to 
contain no information from the measurements, even from those collected prior to the epoch time. 

Substituting equation (31) into equation (36) now gives^ 


N 


N 


i= 1 


p& = P- KjH, 

Pit = S, (y^^kA s ; t , 


\i=l 

N 


N 


N 


P t* =S*Y K ' E [ u di u dj] K JSt - Y K ' E [UdiWd (t* , to)] In - Sit Y 


i,j = 1 


i= 1 


3 = 1 


N 


N 


- Un - Sit Y K ' H; r* E E Yi(P,to)Y] KJSJ 

V i=i / j = t 

N N N 

5* E i^J K Jst - s,Y K ^ E i u ^<(t*,to)} - E E M^^oK] kjs: 

hj = 1 i=1 3 = 1 


T 

* 1 


(37a) 

(37b) 


(37c) 


§ Correlations between e* and are absent in Refs. 2 and 3 because these references implicitly assumed the epoch time to be 
prior to all the measurement times. 

'We use the notation .(•) = 
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where 

= E [udiUjj] + E [u di w d (u, t 0 )] Hj* + Hi,* E [w rf (f*, i 0 )u^-] . (38) 


We define the process noise covariance integrals 


Q*(ti,tj)= / <h(t*,T)Q(T)<£> T (f*,T)dT, 


Q*j — Q*{t*,tj'), 


(39a) 

(39b) 


the second of which, Q*,j, is the usual equivalent discrete-time process noise covariance for the transition 
from time tj to t*. In evaluating these and many subsequent integrals of process noise, care must be taken 
to observe sign changes from reversing the direction of the integration, especially of integrals involving the 
Dirac delta function. In particular, we note that Q*(U,tj) is positive definite if ti > t j and negative definite 
if ti < tj, assuming that the power spectral density Q(r) is positive definite. Using these definitions and the 
relation <f>(ij,r) = $(ij, we have 


E [u di u dj ] = HiE [w d (ti,U)wJ{tj,U)} H' 


3 


= Hi / $(ti,T)Q(T)6(T - T , )$ T (f i ,T , )dTdT , .£/J 

Jt* Jt* 

Hi,* Q*,min(i,j) ^j,* t* ^ ti,t* < tj 
^ Hi,* Q* )m ax(i,j) H j,* ti ^ t*,tj < t *, 

0 otherwise. 


(40) 


We also have 


E[u di wJ(U,t 0 )] = HiE[w d (ti,t*)wJ(t*,t 0 )} 


= -Hi,* / r )Q(t )6(t — r / )4 >T (U, r^drdr' 

J to J ti 


= < 


0 


^ ti 
ti U . 


(41) 


This expectation is simpler than E 


UdiUjj 


because to is less than t * and ti by assumption. We can combine 
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these to compute T ^ : 


r a = < 


Hi,*Q * ,min(i,j 

) H l* 

for p < p, p < tj, 

Hi,* (Q*,max( 

i,j) — Q*,i — Q*,j ) Hjx ^ 

for p < p, tj < p, 


— Hi,*Q *,min(i,j)Hj^. 

for p > p, t,j < p, 

U,..Q..iHJ. 


for p < p, p > p. 


( 42 ) 


Remarkably, Ty takes the same form for any ordering of p, tj, and p. Putting this together gives 


N fc* 

p+* = -S, J2 + J2 (sjuH^Q^i + Q^H^KjS^ ) , (43) 

hi=i 


i = 1 


where k* is the index of the latest measurement time earlier than p, i.e. p,._ < p < th,+\. If p < min (U,tj), 
then <3*, m in(ij) is negative. The process noise component of the net covariance is always positive, however, 
because of the process noise component in P f + implied by equation (35). Equation (43) reduces to the form 
in Refs. 1-3 when p < p; the second sum is identically zero, and all of the Q* iin in(ij) terms in the first sum 
give positive contributions to P+* . 


The state estimation errors propagate to any other time according to equation (12). Therefore, the various 
components of the error covariances at such other times arc given by 1 1 


P a (t) = <&(PP)P+4> T (PP) (44a) 

P v (t) = ^(t,QP+^ T (t,Q (44b) 

P w (t) = $(pp)P+,$ T (pp) + E[w d (pp)wJ(pp)] 

+ $(pp) E [(e+» + e+)wj(t,p)] + {$(pp) E [(e+* + e+ )wj(pp)] } T . (44c) 


We find that 


E[w d (pp)wJ(pp)] 


<b(p p)Q*(max(p p), min(p p))4> T (p p), 


( 45 ) 


I Remember that contains the correlations with the process noise in the a priori error at time t*. 
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and, with equation (35), 


N 


N 


E [( e E + e i) w d (i, Q] = In-S*^2 KiH, i t * E [w d (2*, 2 0 )wJ (2, 2*)] - S* ^ FQE [u dl wj (2, 2*)] 


2=1 


2=1 


= < 


-Q*(2*,2) + 5* Y,iLk+ KiHi i# Q*(imn(U,ti),t) 4> T (2,2*) 2 < 2*, 


(46) 


~S* Eil fe+ ifi-ffi,*Q*(min(2, 2j), 2*)4> T (2, 2* 


2 ^ 2 * , 


where A: + is the index of the earliest measurement time later than either 2 or 2* , i.e. 2m _ i < min (2* , 2) < 2m . 
This gives 

P w {t ) = $(2,2*)[5+* “ ) + 5(2,2*) + 5 T (2,2*)]$ T (2,2*), (47) 


where 


N 

5(2,2*) = <?* ^ 2Tj5j i *Q*(min(2*, 2i), min(2, 2j)). 

2=/v+ 


Note that Q* (2* , 2) is positive if 2 < 2* and negative if 2 > 2* . 


(48) 


Formal Covariance and Delta Covariances 

The estimation system does not take the consider states into account, possibly out of ignorance but more 
likely by deliberate choice to simplify its structure, so it updates an n s -dimensional solve-for state error 
vector. Denoting the a priori estimate of this state error vector at the epoch time by e*i 0 and the updated 
value incorporating all N measurements by e^ N , we have in parallel with equation (30) 

N ( N \ N 

e *| n = e*| 0 - ^2 Ki (5j j *e*| 0 + v, : + u* : ) = I I n ~ ^ KiP,,* e*| 0 - ^ K t {v, + u di ). (49) 

i = 1 V i= 1 J i= 1 

This equation lacks the projection operator onto the solve-for space, because the estimator only deals with 
the solve-for parameters. In this equation 


Hi^ = (50) 

where //, is the in, x n s measurement sensitivity matrix with respect to the solve-for states, <f>(2, ,2*) an 
n s x n s transition matrix of the solve-for states, v, is a Gaussian white measurement noise, and u di is 
defined analogously to equation (29). These arc all the imprecise values modeled by the estimator, and arc 
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not assumed to be equal to the projection onto the solve-for space of their true values.** Of course, an 
estimator that more closely approximates the truth can be expected to provide better estimates. The same 
gain matrices appeal - in the equations for the true and formal covariances, however, because the true gain 
matrix is the one provided by the estimator. 


In addition to providing a solve-for parameter estimate, the estimator will also compute an n s x n s 
covariance matrix estimate, known as the formal error covariance. 


P*\N = E 


e *|lV e *|A r 


= P+ + P+ + P+ 

M a* 1 M v* 1 w* i 


(51) 


Following a path parallel to that leading to equations (37) and (43) gives 


N 


N 


i= 1 


Pa* = In. - E K > H > * ) P *|0 ( ^ - E K ' H >- 

V i = i 

N 

p+ = E^w 


i = 1 


N 


p- = - 


E KiHi,*Q*,min(i,j)Hj,*Kj + E + Q*,jHf*K 7 


(52a) 

(52b) 

(52c) 


i,j= 1 


i= 1 


In these equations f( | 0 is the a priori covariance, II, = E [vjvT] is the measurement error covariance, 


Q*{U,tj) = / 4>(f*,r)Q('r)^ >T (^,'r)dT, 
Q*,j —— Q*(t*,tj) 


(53a) 

(53b) 


with Q(t ) an n s x n s power spectral density of the process noise assumed by the estimator. As was noted 
above, none of these quantities is assumed to be equal to the projection onto the solve-for space of its 
corresponding true value. 


**In the same spirit, 4 > ss (ti, t*) and H a i in Ref. 1 can be replaced by 4>(ti, t*) and fT,*, respectively. 
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The formal covariance is propagated to other times by the analogs of equations (44a), (44b), and (47) 


P(t) — P a (t) + P v (t ) + P w (t) (54a) 

P a (t) = $(t,QP+$ T (t,U) (54b) 

P v (t) =$(t,U)P+$ T (t,Q (54c) 

Pw(t) = - Q*(t*,t) + B(t,t *) + i? T (M*)]<l T (t,t*), (54d) 

where 

N 

B(t, t*) = ^2 (55) 

i=k+ 


The full covariance matrix is needed for computation, but in most cases we are only interested in its 
solve-for part 

Pssit) = E Mf)eJ(i)] = SP a (t)S T + SP v (t)S T + SP w (t)S T . (56) 

This will generally differ from the formal covariance due to estimator mistuning and/or neglect of consider 
parameters in the formal covariance. We express the difference as 

Pssit) = P(t) + A P a (t) + A P v (t) + A P w (t), (57) 

where the delta covariance matrices are 

A P a = SP a S T -P a , A P V = SP V S T -P V , and A P w = SP W S T - P w (58) 

The formal covariance matrices can provide either underestimates or overestimates of the errors, so the delta 
covariance matrices can be positive or negative. 


COMPUTATION OF THE OPTIMAL GAIN MATRIX 

The optimal update gain matrix is chosen to minimize the post-update covariance at the epoch time. The 
total covariance can be written, from equations (51) and (52) as 


N N 

P*\n = P* 10 + E KiWijKj - E [ K AAP* 10 - 4m) + (P* 10 - 

i,j = 1 i= 1 


(59) 
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where 


Wij = HiAP*\0 - 4,min(iJ))^> + RiSij = w; 


The forward process noise covariance is defined as 


QZ = 


Q*,i ti ^ i ■■ 


0 if fy . 


and ();, is the Kronecker delta 


1 for i = j 
0 for i / j. 


We now define the symmetric matrix Wn as an N x N matrix of m, x rrij blocks Wif 


W N = 


W n W l2 ... W 1N 
W21 W 22 ■ ■ ■ W2N 


W m Wn 2 ■ ■ ■ WiSTN 


The inverse Wpf can also be written an N x N matrix of rn, x rrij blocks. 


[W^]n [W^] 12 

_! _ fcV fc*] 22 






fcVl [^V2 ... [W^] NN 


and the relations WnWj y 1 = H r v 1 Wn = I imply that 


^WijlW^jk = Y,^N 1 ]i j W j k = S ik i mi . 
3 = 1 i=i 
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Now equation (51) can be written by completing the square as 


N 


P*\N = P * |0 - E HjAP* 10 - Q7,j) 


i,j = 1 


N 


+ E 


N 


Ki ~ E(^*|o - 


i,j = 1 L k= 1 


Wi 




AT 


^ - E( p *io - 


fc=l 


( 66 ) 


If A and B arc real symmetric matrices, we say that A > B if the matrix A — B is positive semidefinite. 
It is clear that 

N 

^ ~ ' ' ' (67) 


P*\n > P * |0 - E ^ - QPpPl* HUP* 10 - 0 

*J =1 

in this sense, with equality if the measurement gain matrix is chosen to be 


*, 3 > 


N 


Ki = E(A|0 - Qtk)HUW^] ki . 


(68) 


k=\ 


We will use this gain matrix to achieve the minimum P* iat, which can be expressed more compactly as 


N 


P*\N = P* |0 - E K 3 Pj,*{P* |0 - Q 

3 = 1 


*, 3 > 


(69) 


Batch Least-Squares Case 

We want to show that our formulation reduces to the usual batch least-squares estimator if the assumed 
process noise is zero. In that case 

W N = P* P*|o Pj + R, (70) 


where 


P* 


Pi,* 
Hn ,* 


(71) 


and P is block-diagonal with the matrices P, along its main diagonal. The Sherman-Morrison-Woodbury 
matrix inversion lemma 13 gives 


ic- 1 = P" 1 - P^P* 



+ PJP-'P*) 




(72) 
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or 


N 


-1 


lWn 1 } k i = Rk 1 6 ki -Rk 1 Hk,* p Jo+Y, H 7* R J ] H ^* H l* R i l ( 73 ) 


3 = 1 


Substituting into equation (68) with QZk = 0 gives 



N 

/ \ 

Hi ~ A|0 

I ns ~^ R l* R k lR k,* 



k = 1 

V 7=1 J J 


HZAj 1 


(74) 


Now writing I n , = (P* |0 1 + Efc=i ^ 1 H k,*) + Ej=i P J.* P j 1 and performing ap- 

propriate cancellations gives 


-l 


N 


K ‘=\ p :,i+T, H kri'k 

7=1 






(75) 


which is just the usual result obtained by minimizing a quadratic loss function with weighted contributions 
from the measurement innovations and the a priori estimate, as in Refs. 1-3. 


Sequential Implementation 

Inverting the very large matrix Wn is a formidable problem in the general case, and the Sherman- 
Morrison-Woodbury lemma is not applicable if process noise is present. The inversion of this large matrix 
can be avoided even when process noise is present by implementing a sequential estimation procedure. In 
developing a sequential estimator, it will be convenient to denote the gain matrices by K t \ k , where the index 
k is the number of measurements included in the update. The estimation errors after processing k — 1 and k 
measurements are given by equation (49) as 

/ k - 1 \ fc-1 

1 ( In s ^ ^ I e*| 0 ^ ^ T U-di); ("76) 

V i= 1 / i = 1 

( k \ k 

In s ~ 'Z y Hj\ k P i,* J eqo — ^ ^ 4" (*</<)• ("77) 

i=l / i=l 

For a sequential method, we need to be able to write 

d *| k (7n s I^k\kHk,*)^*\k— 1 ^k\kfok T U dk) ■ ("78) 
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Substituting equation (76) into equation (78) and comparing with equation (77) gives the consistency re- 
quirement 

K i\k = (4i s “ K k\kHk,*)Ki\k-i i < k. (79) 

It is not difficult to see that complete agreement among all the terms of equations (76)— (78) follows if the 
gain matrices satisfy this condition. It is shown in the Appendix that equation (79) is always satisfied if the 
estimator ignores process noise, but it can only be satisfied in the presence of process noise if QH = 
for i < k, which means that t- t < t*. We will thus assume that A is in the future of all of the measurements, 
which means that we can remove the arrow from every appearance of QA in the following. 


To make further progress, let W k denote the matrix defined by equation (63), but only including k mea- 
surements. We write this in partitioned form as 


Wk 


Wfc_l 



W kk 


( 80 ) 


where 


V k = 


W lk 

W 2k 


W( k -i)k 


(81) 


Letting P * \ k denote the formal covariance incorporating k measurements, equations (60) and (67) (with 
equality) can be used to show that 


o ~ t\k-i)Hl* 


k—l 


E w ki \w^wj k = VfW^Vk, 

i,j= 1 


(82) 


from which it follows that 


Wkk - VfW^Vk = HkAP*\k-i - Q*,k) Hi* + Rk- 


(83) 
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The standard equation for the inverse of a partitioned matrix gives 




(Wk-1 - VkWj^V?)- 1 


Uk 


uz 


(W kk - VfW^Vk)- 1 


(84) 


where 


Uk = -W^VkiWkk - VfW^Vk)- 1 = -W^VkiW^] 


ttjt— 1 t r. 1-1 


T-l 


T - ll 


kk- 


(85) 


These relations give the gain matrix K k \ k , which is the only one needed for the sequential estimator, as 


fc-i 


K k \k = J2( p * 10 - Q*,i) HAW k \ k = J2( p * 10 - Q*,i) HUW k \ k + (Pqo - Q*,k) H% JW*" 1 ] 


kk 


i — 1 


i= 1 


= (P*| fc _! - Q*,fc) Hi Mk% k = {P*\k-1 - Q*,k) Hi, H k ^P , | fc _i - Q*, fc ) iZ£* + 74 


i-i 


( 86 ) 


The third equality uses, for i < k, 

k— 1 fc— 1 

fe'iifc = = - E V = E 

j = 1 i =1 

fc-i 

= - E - Q*d) V (87) 

j=i 

and then 

fe-i 

E (A|0 - Q*,i) HUW k \]ijHjA p *\ 0 ~ Q*j) = p *\o - P *\k-1, (88) 


and the fourth equality uses equation (83). Equation (86) expresses K k \ k in terms of only P*\ k - \ and 
quantities available at the measurement time t k . Completing the development of the recursive estimator 
requires a recursion relation for P, i k . Using the consistency condition on the gain matrix, equation (79), in 
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equation (69) gives 


k k— 1 

P*\k A|0 — ^ ' ^j\k HjAP*\ 0 — Q*,j) P*]0 — ^k\k H k ,*{P * |0 — Q*,k) — ^ ^ P- j\k 4>(A|0 — Q*j) 

3 = 1 J=1 

H - H k ,*Q*,k 

(2n s P^k\k Pk,*)P*\k— 1 “6 Kk\k H k ,*Q*,k 1 Hk,* (P*\k — 1 Q*,k)- (89) 

It is not surprising that the final form strongly resembles equation (69) with an update from only the last 
measurement and with P*\ k -i as the a priori covariance. 

Relation to the Kalman Filter 

If we define 

P*,k = P*\k-i ~ Q*,ki (90) 

then equation (86) becomes 

K k \ k = P*, fc Hi, (P M P*, fc Hi, + R k )~ 1 (91) 

and equation (89) becomes 

P*,k + 1 — (-fn s H-k\ k Hk,*)P*,k 4“ Q*,k Q*,k + 1 — ( P , Pfc | k Hk,*)P*,k "h Q* (P+1 1 tk)> (92) 

which are equivalent to the equations of the epoch-state estimator used by Montenbruck. 14 It should be 
pointed out, though, that P, k is not the covariance of any obvious error vector. 

If we further make the identification 

P*,k = ^(t*,tk)Pk\k-i^ T (t*,tk) (93) 

and use 

(94) 


fc-i 

= (Jn a — R-k\k H k ,*) P*\o — "y ^ Hj\k— i HjAP*\o — Q*,j) 

3 = 1 
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then equations (91) and (92) become 


(95) 


K k\k = &{t*,tk)Kk 

and 

P k +l\k — ^K^fe+1 j tk) (-^n s K k Hk)P k \k—\^ T ^((fc+l) ^*)Q*(^fc+l) (7fc+l) (*)j (96) 

where 

K k = P k \ k -^{H k Pk \k-i HI + 4)" 1 - (97) 

Notice that H /,. appears in these equations rather than H k ,*- Equation (96) can be written as the usual 
two-step Kalman filter process 

p k \k = ( Ins - Kk H k )P k \k-l (98) 

x „ rtk + 1 „ 

-Pfc+i|fc = $(tk+i,tk)Pk\k® T (tk+i,tk) + / $(tfc+i) "7" )Q(r) < l )T (tfc+i, r)dr, (99) 

Jtk 

where the expression for the process noise follows from equation (53). It is now easy to recognize P k \h-i 
and P k i k as the pre- and post-update covariances of the estimation error at the measurement time t, k . 

SUMMARY AND CONCLUSION 

We have updated and extended previous work on consider covariance analysis to improve the treatment 
of process noise when the update time is not prior to all of the measurements. Consider covariance analysis 
treats estimators that account for only a subset of the parameters, possibly out of ignorance but more likely 
by deliberate choice to simplify their structure. The covariance analysis, however, considers the effect on the 
covariance of all the variables, including those ignored by the estimator. Our analysis explicitly partitions 
the full covariance and solve-for covariance into three subspaces containing the influence of measurement 
noise, process noise, and a priori errors, respectively. Explicitly exhibiting the a priori component of the 
covariance is useful because a priori information about observable parameters is “forgotten” as more data 
arc processed, providing a useful test of observability. 

We have found an expression for the gain of an optimal epoch-state estimator incorporating a model of 
process noise, and showed that this reduces to the usual batch estimator if process noise is ignored when 
computing the measurement update gains. Our formulation agrees with our previous analysis of batch 
filtering methods if the epoch time is prior to all measurement times, but differs if the epoch time is amid 
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the data or in the future of all data. If the epoch time is no earlier than the time of the last measurement 
processed, a sequential implementation of the filter is possible, and we show that this is basically equivalent 
to an extended Kalman filter. 

APPENDIX: CONSISTENCY CONDITION FOR SEQUENTIAL IMPLEMENTATION 

We will first show that the consistency condition requires that either the estimator ignores process noise 
or that ti < t* for i < k. We will use a very simple model with a one-component state vector and two scalar 
measurements with iT,,* = H 2 ,* = 1. The blocks 


mj = (A 10 


)) T A, 8 % j 


(100) 


of Wk are scalars, so we have 


W-, 


-1 


1 

W 11 



1 

W n W 2 2 - W[ 2 


W 22 

-W 12 


-W 12 

W n 


The gains are given by equation (68) as 


-At 1 2 
A 212 


A|o ~ Q 7 a 

W 11 

( A |0 - Q ^) W 2 2 - ( A |0 - Q ^ 2 )^12 
W n W 2 2 - Wl 2 

(A|o - Qt-iWu - (A | 0 - Qt„ i)Wi 2 

U , ,1122 \V? 2 


The consistency condition, equation (79) then gives 


(101a) 

(101b) 


(102a) 

(102b) 

(102c) 


(A|o - Q^i)[WnW22 ~ Wl ■ - (A|o - Q??)W: n + (A|o - Q£i)Wu\ = 

Wn[(A|o - Q 7 ,i)w 2 2 - (A | 0 - Q 7 , 2 ) w 12] (103) 

Canceling the common term of (P*|o — \ ) H i 1 W 22 front the two sides of the equation and rearranging 

gives 

[(A|0 - QJ 2 )W 11 - (A|0 - Qm)Wi 2 ](A|0 - QA - W 12 ) = 0 (104) 
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The quantity in square brackets is certainly not zero, because there is nothing to cancel the term R\ in W\\. 
Thus we must have 

o = P*|0 - Q7,1 - W 12 = 4,1 - Q^i (105) 

This equation is satisfied if either the estimator ignores process noise entirely or if t\ < t t . In the latter case 
it follows by by induction that we must have < f* for all i < k. 


The next task is to show that the consistency condition holds if = Q*,i f° r i < k. We will show this in 
the general case. Equation (68) gives, with the subscripting notation used for the sequential implementation, 


k^i\k 


k k— 1 

£(A|o -Q*,j)HjA^k% = l^( p *\o-Q*j)HUWk% + ( P *\o-^ (106) 

3 = 1 i =1 


We use equation (84) and the Sherman-Morrison-Woodbury matrix inversion lemma write the upper left- 
hand corner of Wj~ 1 in the form 


(W k -1 - VkW^V?)- 1 = w-\ + W^\V k (W kk - VfW^Vk)- 1 ^ 

= + W k \V k [W^jkkV? w£ lt 


(107) 


which is made up of the blocks 


l w k % = I'Wib. + [w^v^ \w ^ ] u k t ITT-ik 


(108) 


We also have 


\W k l ]ki = [m = -[W^}kk [Vk T Ufcji 


— ll . . fT/T tjt i 


(109) 


With these substitutions. 


k — 1 k — 1 

Ki\k = J2( p * 10 - Q*j)HJAWk-ihi + E^io - Q^HjAw^Vk], [w-\ k [v? w^]i 

3 = 1 i =1 


- (P H0 - Q*,k)H k *\W^ l ]kk [Vk W^]i. 


(HO) 


The first term on the right side of this equation is easily recognized to be K t \ k - \ ■ The sum over j in the 
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second term can be written as 

k — 1 k — 1 

E( p * 10 - Q*j)H*[Wi\v k \j = Y (A|0 - 

J=1 j,Tn=l 

k - 1 

= E (A|0 - 

j,m=l 

k- 1 

= £ (A|0 - [W^ljmHm,* (P*|0 - = (A|0 ~ (HD 

j,m=l 

where equation (67) (with equality) gives the last line. Putting all this together gives 

tfii* = Ki\k -1 + (A |0 - [^ T WT-Ji - (A |0 - Q*,k)HlM^] kk [v k T w k \]i 

k - 1 

= ifiifc-t - (A| fc -t - 

3 = 1 

fc-1 

= #*!*_! - (P*|fc-i - Q*,k)HlAW k % k Y W kj [W^]ji 

3 = 1 

fc-1 

= ^Ife-r - (A|fc-1 - 

i= 1 

= [-^n s — (-f*|fc— 1 — Q*,k)^k,* \^k ]fcfc-^fc,*]-^Q|fc— 1> = (^n s — -^fc|fc-^fc,*)-^i|fc— 1 (112) 

recognizing equation (86) in the last line. This is the desired relation. 
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